## Estimation of a Multinomial Logit Model by Maximum Likelihood using mnlogit and Fixed point subroutine with continuous weightingfunction
## by Jose-Alberto Guerra
##    Created   04/2014
##    Modified  17/06/2019
##    Uses: R version 3.1.2 (install old version of R, by holding down the Control key during the launch of RStudio you can use 64-bit 3.1.2 version)
#PooledNAPPIDnewwcontinuous.RData comes from Mlogit_NAPPHeterov4
#NAPPXlistnewwcontinuous.RData comes Mlogit_NAPPHeterov4
#coordslistNAPPnewwcontinuous.RData comes Mlogit_NAPPHeterov4
#to deal with spatial data visualitaion in R see http://spatial.ly/r/
## Where we guarantee that KS probabilities are within [0,1]
## We guarantee inner loop (fixed point expectations) convergence using tol 10e-10 for the sup-norm stopping criterion
## We follow KS 2012 and KS 2018 using outer loope iterations k=50 and guaranteeing inner loop convergence (see above)
###########################################################################
rm(list=ls(all=TRUE))
#########################
####    Preamble     ####
#########################
##Functions
#codes <- "C:/Users/ja.guerra/Dropbox/NetworkParish2016/Occupational_Choice/Codes/RCodes/EEACodes/"
codes <- "C:/Users/ja.guerra/Dropbox/NetworkParish2016/Occupational_Choice/Submission/Final_ReStat/Codes/C_ANALYSIS/" #codes <- "C:/HPC/CODE/"
source(paste(codes,"Mlogit_NAPPHeteroFunctionsv2.R",sep=""))
###########################################################################
#1. DATA
###########################################################################
###########################################################################

set.seed(12345)
#define WD and DATASET
##If normal data set
## root <- "C:/Users/ja.guerra/Dropbox/NetworkParish2016/Occupational_Choice/" #root <- "C:/HPC/"
root <- "C:/Users/ja.guerra/Dropbox/NetworkParish2016/Occupational_Choice/Submission/Final_ReStat/Codes/" #root <- "C:/HPC/"
mapdata <- paste(root,"MOLA/",sep="")
dataOGR <- paste(root,"MOLA/",sep="")
#wddata<-c(paste(root,"Results/EEA/",sep=""),paste(root,"NAPP/Proof2RCoordnew.dta",sep="")) #If Placebo dataset #wddata<-c("C:/HPC/Results/Placebo","C:/HPC/NAPP_address/Proof2RPlacebo.dta")
wddata<-c(paste(root,"Results/ReStat/",sep=""),paste(root,"NAPP/Proof2RCoordnew.dta",sep="")) #If Placebo dataset #wddata<-c("C:/HPC/Results/Placebo","C:/HPC/NAPP_address/Proof2RPlacebo.dta")
setwd(wddata[1])

sink("ExplorenewWContinuous.txt")

#Calling Data see Mlogit_NAPPHeterov3.R
load(file=paste(root,"NAPP/ReStat/PooledNAPPIDnewwcontinuous.RData",sep="" )) 
load(file=paste(root,"NAPP/ReStat/NAPPXlistnewwcontinuous.RData",sep="" )) 
load(file=paste(root,"NAPP/ReStat/coordslistNAPPnewwcontinuous.RData",sep="" )) 

#### From here onwards it follows Mlogit_NAPPHeterov4.R from line 339
coordslist <- coordslistNAPP


##############################################
#1. Robustness weighting matrices 
## W continuous distance
## Per paris, list of neighbours
dallneighbours <- lapply(coordslistNAPP, withinAllnei) 

##Per parish, list of neighbours row standardised (this is what we use in our Mlogit_NAPPHeterov3)
wz <- lapply(dallneighbours, function(z) nb2listw(z, glist=NULL, style="W", zero.policy=TRUE)) 
#wz[[1]],...,wz[[G]]
# Notice wz[[g]]$weights[[1]],...,Notice wz[[g]]$weights[[n_g]]

##Per parish, list of distances between neighbours
parishdistances <- function(x,y){
  nbdists(x,y,longlat=NULL)
}
dz <- mapply(parishdistances,dallneighbours,coordslistNAPP)
#dz[[1]],...,dz[[G]]
#Notice dz[[g]][[1]],...,dz[[g]][[n_g]]
#############
### some descriptive statistics
d50neighbours <- lapply(coordslistNAPP, within50nei )
descdz <- c(1:nrow(PooledNAPPID)) 
descn50 <- c(1:nrow(PooledNAPPID)) 
lparish <- 0
for (pp in c(1:length(wz))){
  for(ii in c(1:length(wz[[pp]]$weights))){
    descdz[lparish+ii] <- mean(dz[[pp]][[ii]])
    descn50[lparish+ii] <- length(d50neighbours[[pp]][[ii]])
  }
  lparish <- lparish+length(wz[[pp]]$weights)
}
#Distance between pairs
sapply(list(mean,median,sd,min,max),function(fun,x) fun(descdz))
#Number of neighbors within 50 meters
sapply(list(mean,median,sd,min,max),function(fun,x) fun(descn50))
#########
##Per parish, list of neighbours row standardised (this is what we use in our Mlogit_NAPPHeterov3)
walllist <- wz
for (pp in c(1:length(walllist))) {
  #Replacing count weights by exp(-distance weights)/sum(exp(-distance weights)), row normalization
  walllist[[pp]]$weights <- lapply(dz[[pp]],function(x) exp(-x)/sum(exp(-x)) )  
  }


dneighbours <- dallneighbours
wlist <- walllist

##############################################

#Combining lists
NAPPXwlist <- mapply(weighting1,NAPPXlist,wlist) 
NAPPXlist <- lapply(NAPPXlist,function(x) x <- as.data.frame(x)) #to convert  matrix to data.frame
pooledNAPPXlist <- rbindlist(NAPPXlist)
pooledNAPPXlist <- cbind(PooledNAPPID, pooledNAPPXlist)

NAPPXwlist <- lapply(NAPPXwlist,function(x) x <- as.data.frame(x)) #to convert  matrix to data.frame
pooledNAPPXwlist <- rbindlist(NAPPXwlist)
attr(pooledNAPPXwlist,"class") <- "data.frame"


pooleddatalist <- cbind(pooledNAPPXlist[ !(colnames(pooledNAPPXlist) %in% paste("class",c(min(pooledNAPPXlist$choice):max(pooledNAPPXlist$choice)),sep="")) ],pooledNAPPXwlist)
colnames(pooleddatalist)[colnames(pooleddatalist) %in% paste("class",c(min(pooledNAPPXlist$choice):max(pooledNAPPXlist$choice)),"w",sep="") ] <- paste("sw.",c(min(pooledNAPPXlist$choice):max(pooledNAPPXlist$choice)),sep="")
datalist <- split(pooleddatalist, list(pooleddatalist$IDSocial))

#Poduce the long format and pooled data
datalistlong <- lapply(datalist,function(x) tolong(x,"choice",c("sw.")))
pooleddata <- rbindlist(datalistlong)

#Produce long data for homogenoeus beliefs
pooleddatalisthomo <- pooleddatalist
pooleddatalisthomo[grep("w", names(pooleddatalisthomo),value=TRUE)] <-  ddply(pooledNAPPXlist, .(IDSocial), sapply, FUN = 'witaver')[c("SEX","AGE","married",  "nchild",   "nservant", "stayerp",  "stayerc",   "class0",   "class1",   "class2", "class3", "class4", "class5")]

pooleddatahomo <-  tolong(pooleddatalisthomo,"choice",c("sw."))

#To save data paste(root,"Mola Data/",sep="")
save(pooleddata, file=paste(root,"NAPP/ReStat/pooleddataNAPPnewwcontinuous.RData",sep="" )) #Data per parish in long format (sw) pooled together
save(pooleddatahomo, file=paste(root,"NAPP/ReStat/pooleddataNAPPhomonewwcontinuous.RData",sep="" )) #Data per parish in long format (sw) pooled together for homogeneous beliefs
save(pooleddatalist, file=paste(root,"NAPP/ReStat/pooleddatalistNAPPnewwcontinuous.RData",sep="" )) #Data per parish in wide format (sw.1,sw.2,...) pooled together
save(pooleddatalisthomo, file=paste(root,"NAPP/ReStat/pooleddatalistNAPPhomonewwcontinuous.RData",sep="" )) #Data per parish in wide format (sw.1,sw.2,...) pooled together for the homogeneous case
save(datalist, file=paste(root,"NAPP/ReStat/datalistNAPPnewwcontinuous.RData",sep="" )) #Lists of data per parish in wide format (sw.1,sw.2,...)
#save(coordslist, file=paste(root,"NAPP/ReStat/coordslistNAPPnewwcontinuous.RData",sep="" )) #Lists of coordinates per parish
save(dneighbours, file=paste(root,"NAPP/ReStat/dneighboursNAPPnewwcontinuous.RData",sep="")) #Lists of neigbours identities per parish
save(wlist, file=paste(root,"NAPP/ReStat/wlistNAPPnewwcontinuous.RData",sep="")) #Sparse Lists of neighbours per parish


sink()

#########################################################################
###2. ESTIMATION######
########################################################################
#  rm(list=ls(all=TRUE))
#  set.seed(12345)
#  #root <- "S:/NetworkParish/HPC/"
# root <- "C:/Users/ja.guerra/Dropbox/NetworkParish2016/Occupational_Choice/" #root <- "C:/HPC/"
# setwd(paste(root,"Results/ReStat/",sep=""))
# #codes <- "C:/Users/ja.guerra/Dropbox/NetworkParish2016/Occupational_Choice/Codes/RCodes/EEACodes/" #codes <- "C:/HPC/CODE/"
# codes <- "C:/Users/ja.guerra/Dropbox/NetworkParish2016/Occupational_Choice/Submission/ReStat_R&R/Codes/C_ANALYSIS/" #codes <- "C:/HPC/CODE/"
# source(paste(codes,"Mlogit_NAPPHeteroFunctionsv2.R",sep=""))
# #Load Data
# 
# load(file=paste(root,"NAPP/ReStat/pooleddataNAPPnewwcontinuous.RData",sep="" )) # pooleddata, Data per parish in long format (sw) pooled together
# load(file=paste(root,"NAPP/ReStat/pooleddataNAPPhomonewwcontinuous.RData",sep="" )) # pooleddatahomo, Data per parish in long format (sw) pooled together for hmogenous
# load(file=paste(root,"NAPP/ReStat/pooleddatalistNAPPnewwcontinuous.RData",sep="" )) #pooleddatalist, Data per parish in wide format (sw.1,sw.2,...) pooled together
# load(file=paste(root,"NAPP/ReStat/pooleddatalistNAPPhomonewwcontinuous.RData",sep="" )) #pooleddatalisthomo, Data per parish in wide format (sw.1,sw.2,...) pooled together
# load(file=paste(root,"NAPP/ReStat/datalistNAPPnewwcontinuous.RData",sep="" )) #datalist, Lists of data per parish in wide format (sw.1,sw.2,...)
# load(file=paste(root,"NAPP/ReStat/coordslistNAPPnewwcontinuous.RData",sep="" )) #coordslist, Lists of coordinates per parish
# load(file=paste(root,"NAPP/ReStat/dneighboursNAPPnewwcontinuous.RData",sep="")) #dneighbours, Lists of neigbours identities per parish
# load(file=paste(root,"NAPP/ReStat/wlistNAPPnewwcontinuous.RData",sep="")) #wlist, Sparse Lists of neighbours per parish
# 
# ###########################
# ###Homogeneous Beliefs (i.e. group interactions, no correlated effects and fixed point beliefs subroutine)
# ###########################
# AM <- 0 #if using Aguirregabiria and Mira==1, 0 if using Kasahara and Shimotsu
# 
# #if( AM==1 ){  rootsave <- paste(root,"Results/EEA/",sep="") } else   {rootsave <- paste(root,"Results/EEA/KS/",sep="")}
# if( AM==1 ){  rootsave <- paste(root,"Results/ReStat/",sep="") } else   {rootsave <- paste(root,"Results/ReStat/KS/",sep="")}
# sink(paste(rootsave,"NAPPEstimationHomonewwcontinuous.txt",sep=""))
# print("Homogenous beliefs: group interations, no correlated effects and fixed point subroutine")
# 
# #define which variables we include in the estimation
# xvars <- c("AGE","married",  "nchild",   "nservant", "stayerp") #, "stayerc" "SEX"
# wvar <- "sw"
# 
# ## Naive Estimation Homogeneous
# ###################
# print("Naive Estimation")
# options(max.print=1000000)
# gc()
# #1. x_i
# fformula <-mFormula(as.formula(paste("choice ~ 1 | ", paste(xvars, collapse= "+"),"|", wvar ,sep="")))
# ml.w <- mnlogit(fformula, pooleddatahomo,"alt")
# save(ml.w,file=paste("Mlwhomowcontinuous.RData",sep=""))
# print("w, x_i")
# print(summary(ml.w))
# #2. x_i, x_p
# ##w
# fformula <-mFormula(as.formula(paste("choice ~ 1 | ", paste(paste(xvars, collapse= "+"),"+",paste(xvars, "w",collapse= "+",sep="")),"|", wvar ,sep="")))
# ml.om.w <- mnlogit(fformula, pooleddatahomo,"alt")
# save(ml.om.w,file=paste("Mlomwhomowcontinuous.RData",sep=""))
# print("w, x_i x_p")
# print(summary(ml.om.w))
# #3. x_i, x_p, t_b
# ##w
# fformula <-mFormula(as.formula(paste("choice ~ 1 | ", paste(paste(xvars, collapse= "+"),"+",paste(xvars, "w",collapse= "+",sep=""),"+","factor(IDCivil)"),"|", wvar ,sep="")))
# ml.m.w <- mnlogit(fformula, pooleddatahomo,"alt")
# save(ml.m.w,file=paste("Mlmwhomowcontinuous.RData",sep=""))
# print("w, x_i x_p t_b")
# print(summary(ml.m.w))
# 
# 
# ## Nested Pseudo likelihood Estimation Homogeneous
# ###################
# 
# #define which variables we include in the estimation and which formula we use
# xvars <- c("AGE","married",  "nchild",   "nservant", "stayerp") #, "stayerc" "SEX",
# wvar <- "sw"
# fformula <-mFormula(as.formula(paste("choice ~ 1 | ", paste(paste(xvars, collapse= "+"),"+",paste(xvars, "w",collapse= "+",sep=""),"+","factor(IDCivil)"),"|", wvar ,sep="")))
# 
# #Outer Iteration: Maximum likelihood
# alter <- sort(unique(pooleddatahomo$alt))
# L <- length(alter)
# it0 <- 1
# #maxiter0 <- 600
# #zero <- 1E-3
# maxiter0 <- 50
# zero0 <- 1E-8
# Jeit <- matrix(rep(0,maxiter0*L),nrow=maxiter0) #Keep track of the coefficients
# ndJeit <- L
# 
# #To keep original pooleddata to implement Weidner (2012) method
# pooleddatahomo0 <- pooleddatahomo
# pooleddatalisthomo0 <- pooleddatalisthomo
# 
# #to define outer iteration of expectations
# ndexpectito <- L
# ndexpectitoall <- matrix(0,nrow=maxiter0) #Keep track of the tolerance level
# expeco <- subset(pooleddatalisthomo,select=paste(wvar,".",alter,sep=""))
# 
# #while (it0<=maxiter0 & ndJeit > zero) {
# while (it0<=maxiter0 & ndexpectito > zero0) {
#   print(paste("ML iteration: ",it0, sep=""))
#   #Estimating
#   if(it0==1) {
#       load("Mlmwhomowcontinuous.RData")
#       ml.m.w0 <- ml.m.w
#     } else  ml.m.w <- mnlogit(fformula, pooleddatahomo, "alt")
# 
#   #Define starting conditions, if we want to give them staring values based on theta
#   #if (it0==1) ml.m.w$coeff[paste(wvar,":",alter,sep="")] <- runif(L,2,4)
# 
#   #Keeping coeff
#   Jeit[it0,] <- summary(ml.m.w)$coeff[paste(wvar,":",alter,sep="")]
#   print(Jeit[it0,])
#   ndexpectitoall[it0,] <- ndexpectito 
#   print(ndexpectitoall[it0,])
#   #Inner iteration, fixed point in beliefs
#   it <- 1
#   zero <- 1E-8
#   #maxiter <- 1 #if iterating a la Aguirregabiria and Mira (2007) or a la Weidner (2012)
#   maxiter <- 100 #If iterating with the fixed point a la Lee Lin Liu (2014)
#   expectit <- matrix(rep(0,maxiter*L),nrow=maxiter)
#   ndexpectit <- L
#   #subsetting sample to
#   expec <- subset(pooleddatalisthomo,select=paste(wvar,".",alter,sep=""))
#   #iteration of beliefs
#   while (it<=maxiter & ndexpectit > zero) {
#     print(paste("Inner Iteration ",it,sep=""))
#     #keepint track of expectations
#     expectit[it,] <- apply(expec,2,mean)
# 
#     alpha <- .1 #Kasahara and Shimotsu (2012)
#     #Predict probabilities on the pooled data list
#     tempprob <- predict(ml.m.w,pooleddatahomo,probability=TRUE)
#     #Update Fixed point
#     for(aa in sort(alter, decreasing=TRUE) ){
#       ap <- which(alter==aa)
#       if (aa != min(alter)) { #to guarantee KS probabilities are within [0,1] bounds
#         if (AM==1) { 
#           pooleddatalisthomo[,paste(wvar,".",aa,sep="")] <- ml.m.w$probabilities[,ap] #Aguirregabiria and Mira (2007)
#         } else {
#           #pooleddatalisthomo[,paste(wvar,".",aa,sep="")] <- ((ml.m.w$probabilities[,ap])^alpha)*((ml.m.w0$probabilities[,ap])^(1-alpha)) #Kasahara and Shimotsu (2012)
#           pooleddatalisthomo[,paste(wvar,".",aa,sep="")] <- ((tempprob[,ap])^alpha)*((ml.m.w0$probabilities[,ap])^(1-alpha)) #Lee and Liu (2014)
#         }
#       } else {
#         if (AM==1) { 
#           pooleddatalisthomo[,paste(wvar,".",aa,sep="")] <- 1-apply(ml.m.w$probabilities[,paste(alter[alter!=min(alter)],sep="")],1,sum) 
#         } else {
#           pooleddatalisthomo[,paste(wvar,".",aa,sep="")] <- 1-apply(((tempprob[,paste(alter[alter!=min(alter)],sep="")])^alpha)*((ml.m.w0$probabilities[,paste(alter[alter!=min(alter)],sep="")])^(1-alpha)),1,sum) #Kasahara and Shimotsu (2012)
#         }
#       }
#       #pooleddatalisthomo[,paste(wvar,".",aa,sep="")] <- predict(ml.m.w,pooleddatahomo,probability=TRUE)[,ap]  #Aguirregabiria and Mira (2007)
#       #Weidner
#       #pooleddatalisthomo[,paste(wvar,".",aa,sep="")] <- predict(ml.m.w,pooleddatahomo0,probability=TRUE)[,ap]  #Weidner
#     }
# 
#     #``spatially'' weighting
#     pooleddatalisthomo[paste(wvar,".",alter,sep="")] <-  ddply(pooleddatalisthomo, .(IDSocial), sapply, FUN = 'witaver')[paste(wvar,".",alter,sep="")]
# 
#     #updating expectations
#     expec1 <- subset(pooleddatalisthomo,select=paste(wvar,".",alter,sep=""))
# 
#     #Updating pooleddata (long format) #Another option pooleddatahomo <-  tolong(pooleddatalisthomo,"choice",c(paste(wvar,".",sep="")))
#     for(aa in alter){
#       pooleddatahomo[pooleddatahomo$alt==aa,][[paste(wvar,sep="")]] <- pooleddatalisthomo[[paste(wvar,".",aa,sep="")]]
#     }
# 
#     if (it>=2) ndexpectit <- max(as.matrix(abs(expec1-expec))) # the maximum of all matrix entries #ndexpectit <- norm(as.matrix(expec1-expec),type="I")  #element by element sum(distr(expec,expec1))
#     print(paste("Max abs diff inner: ",ndexpectit,sep=""))
#     expec <- expec1
#     it <- it+1
#   }
#   if (it<=maxiter & ndexpectit < zero) print(paste("Fixed point in beliefs convergence reached at iter: ",it-1,sep="")) else print(paste("Fixed point in beliefs one step update it: ",it,sep="")) #print(paste("Fixed point in beliefs FAILED to converge after it: ",it,sep=""))
#   
#   #updating expectations
#   expec1o <- subset(pooleddatalisthomo,select=paste(wvar,".",alter,sep=""))
#   
#   #if (it0>=2) ndJeit <- sum(abs(Jeit[it0,]-Jeit[it0-1,]))
#   if (it0>=2) {
#     #ndexpectito <- norm(as.matrix(expeco-expec1o),type="I")
#     ndexpectito <- max(as.matrix(abs(expec1o-expeco))) # the maximum of all matrix entries
#     }#element by element ndexpectito <- sum(distr(expeco,expec1o))
#   print(paste("Diff Abs max ",ndexpectito,sep=""))
#   expeco <- expec1o
#   #To keep the previous probabilities to perform Kasahara and Shimotsu
#   ml.m.w0 <- ml.m.w
#   it0 <- it0+1
# }
# print(paste("Outer ML Homogeneous iter converged at iter: ",it0-1," using weights: ",wvar,sep=""))
# Jeit[it0-1,] <- summary(ml.m.w)$coeff[paste(wvar,":",alter,sep="")]
# print(Jeit[it0-1,])
# ndexpectitoall[it0-1,] <- ndexpectito 
# print(ndexpectitoall[it0-1,])
# #Saving last iteration results
# save(ml.m.w,file=paste(rootsave,"PMLmwhomofinalwcontinuous.RData",sep=""))
# print(summary(ml.m.w))
# 
# #Plot convergence in coefficientes
# #coefficients
# colnames(Jeit)<-c(paste("J.",alter,sep=""))
# dataJeit<-as.data.frame(Jeit[c(1:(it0-1)),])
# dataJeit$iter <- c(1:nrow(dataJeit))
# meltdataJeit<-melt(dataJeit,id="iter")
# cairo_ps(paste(rootsave,"JePMLmwhomofinalwcontinuous.eps",sep=""), width=7, height=7) #to keep transparency from the intervals
# print(ggplot(meltdataJeit,aes(x=iter,y=value,colour=variable,group=variable)) + geom_line()+ ylab(expression(paste(J[y]^{t})))+xlab("t"))
# dev.off()
# save(dataJeit,file=paste(rootsave,"JePMLmwhomofinalwcontinuous.RData",sep=""))
# 
# 
# 
# #Plot convergence in Infinity Norm
# #coefficients
# colnames(ndexpectitoall)<-"AbsMax"
# dataInftyNorm<-as.data.frame(ndexpectitoall[c(3:(it0-1))])
# colnames(dataInftyNorm)<-"AbsMax"
# dataInftyNorm$iter <- c(1:nrow(dataInftyNorm))
# meltdataInftyNorm<-melt(dataInftyNorm,id="iter")
# cairo_ps(paste(rootsave,"InftyNormMLmwhomofinalwcontinuous.eps",sep=""), width=7, height=7) #to keep transparency from the intervals
# print(ggplot(meltdataInftyNorm,aes(x=iter,y=value,colour=variable,group=variable)) + geom_line()+ ylab(expression(InftyNorm))+xlab("t"))
# dev.off()
# save(dataInftyNorm,file=paste(rootsave,"InftyNormMLmwhomofinalwcontinuous.RData",sep=""))
# sink()


###########################
###Heterogenous Beliefs (i.e. ``network'' interactions, correlated effects and fixed point beliefs subroutine)
###########################
rm(list=ls(all=TRUE))
set.seed(12345)
#root <- "S:/NetworkParish/HPC/"
root <- "C:/Users/ja.guerra/Dropbox/NetworkParish2016/Occupational_Choice/Submission/Final_ReStat/Codes/" #root <- "C:/HPC/"
setwd(paste(root,"Results/ReStat/",sep=""))
codes <- "C:/Users/ja.guerra/Dropbox/NetworkParish2016/Occupational_Choice/Submission/Final_ReStat/Codes/C_ANALYSIS/" #codes <- "C:/HPC/CODE/"
source(paste(codes,"Mlogit_NAPPHeteroFunctionsv2.R",sep=""))
AM <- 0 #if using Aguirregabiria and Mira==1, 0 if using Kasahara and Shimotsu
#if( AM==1 ){  rootsave <- paste(root,"Results/EEA/",sep="") } else   {rootsave <- paste(root,"Results/EEA/KS/",sep="")}
if( AM==1 ){  rootsave <- paste(root,"Results/ReStat/",sep="") } else   {rootsave <- paste(root,"Results/ReStat/KS/",sep="")}
#codes <- "C:/Users/ja.guerra/Dropbox/NetworkParish2016/Occupational_Choice/Codes/RCodes/EEACodes/" #codes <- "C:/HPC/CODE/"
#source(paste(codes,"Mlogit_NAPPHeteroFunctions.R",sep=""))

#######
#Load Data

load(file=paste(root,"NAPP/ReStat/pooleddataNAPPnewwcontinuous.RData",sep="" )) # pooleddata, Data per parish in long format (sw) pooled together
load(file=paste(root,"NAPP/ReStat/pooleddataNAPPhomonewwcontinuous.RData",sep="" )) # pooleddatahomo, Data per parish in long format (sw) pooled together for hmogenous
load(file=paste(root,"NAPP/ReStat/pooleddatalistNAPPnewwcontinuous.RData",sep="" )) #pooleddatalist, Data per parish in wide format (sw.1,sw.2,...) pooled together
load(file=paste(root,"NAPP/ReStat/pooleddatalistNAPPhomonewwcontinuous.RData",sep="" )) #pooleddatalisthomo, Data per parish in wide format (sw.1,sw.2,...) pooled together
load(file=paste(root,"NAPP/ReStat/datalistNAPPnewwcontinuous.RData",sep="" )) #datalist, Lists of data per parish in wide format (sw.1,sw.2,...)
load(file=paste(root,"NAPP/ReStat/coordslistNAPPnewwcontinuous.RData",sep="" )) #coordslist, Lists of coordinates per parish
load(file=paste(root,"NAPP/ReStat/dneighboursNAPPnewwcontinuous.RData",sep="")) #dneighbours, Lists of neigbours identities per parish
load(file=paste(root,"NAPP/ReStat/wlistNAPPnewwcontinuous.RData",sep="")) #wlist, Sparse Lists of neighbours per parish

sink(paste(rootsave,"NAPPEstimationHetenewwcontinuous.txt",sep=""))
options(max.print=1000000)
gc()
print("Heterogeneous  beliefs: ``network'' interactions, correlated effects and fixed point beliefs subroutine")

#define which variables we include in the estimation
xvars <- c("AGE","married",  "nchild",   "nservant", "stayerp") #, "stayerc" "SEX",
wvar <- "sw"

## Naive Estimation Heterogeneous
###################
print("Naive Estimation")
gc()
#1. x_i  
fformula <-mFormula(as.formula(paste("choice ~ 1 | ", paste(xvars, collapse= "+"),"|", wvar ,sep="")))          
ml.w <- mnlogit(fformula, pooleddata,"alt")
save(ml.w,file=paste("Mlwnewwcontinuous.RData",sep=""))
print("w, x_i")
print(summary(ml.w))
#2. x_i, x_p
##w
fformula <-mFormula(as.formula(paste("choice ~ 1 | ", paste(paste(xvars, collapse= "+"),"+",paste(xvars, "w",collapse= "+",sep="")),"|", wvar ,sep="")))          
ml.om.w <- mnlogit(fformula, pooleddata,"alt")
save(ml.om.w,file=paste("Mlomwnewwcontinuous.RData",sep=""))
print("w, x_i x_p")
print(summary(ml.om.w))
#3. x_i, x_p, t_b
##w
fformula <-mFormula(as.formula(paste("choice ~ 1 | ", paste(paste(xvars, collapse= "+"),"+",paste(xvars, "w",collapse= "+",sep=""),"+","factor(IDCivil)"),"|", wvar ,sep="")))          
ml.m.w <- mnlogit(fformula, pooleddata,"alt")
save(ml.m.w,file=paste("Mlmwnewwcontinuous.RData",sep=""))
print("w, x_i x_p t_b")
print(summary(ml.m.w))

#4. x_i, x_p, t_s
##w
fformula <-mFormula(as.formula(paste("choice ~ 1 | ", paste(paste(xvars, collapse= "+"),"+",paste(xvars, "w",collapse= "+",sep=""),"+","factor(IDSocial)"),"|", wvar ,sep=""))) 
ml.mt.w <- mnlogit(fformula, pooleddata,"alt")
save(ml.mt.w,file=paste("Mlmtwnewwcontinuous.RData",sep=""))
#0 When robustness we have to change the ending of the file to accomodate one of the following
# #wexernei COULD BE atleast1 w100list nmovers age1530 w50list atmost10
#1. end of robust
print("w, x_i x_p t_p")
print(summary(ml.mt.w))
sink()

## Nested Pseudo likelihood Estimation Heterogeneous
###################

####Mak sure that wexernei is wcontinuous
#source(paste(codes,"Mlogit_NAPPHeterofformula1Robustv2.R",sep=""))